## Loading required package: knitr
## Loading project configuration
## Autoloading helper functions
##  Running helper script: funRel.R
## Autoloading cache
## Autoloading data
##  Loading data set: DARTallIdentifiers
##  Loading data set: qslAllLarvaInfo
##  Loading data set: qslMetaLarvAndAdultsUnion
##  Loading data set: qslMPeeliiForRelated
##  Loading data set: Report.DMac15.1861
## Munging data
##  Running preprocessing script: 01MungeGeneticsDatacontaminationfix.R
## Loading required package: dplyr
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
##  Running preprocessing script: 02MungeLarvaSnpsByYears.R
##  Running preprocessing script: 03MungeLarvaSnpsByYearCombos.R
a<-suppressPackageStartupMessages({
library(dplyr)
library(igraph)
library(ggplot2)
})

ptm <- proc.time()

This uses iGraph to plot relatedness between larvae within each year so as to identify full sibling pairs(FS) and perhaps half sibling pairs (HS). Unrelated (US) is also calculated by the ‘r’ package ‘related’. However we are interested only in FS at this stage to identify common parents and to assist with determining the distance of larval dispersal.

The same is finally ploted for all three years combined.

First choose state which estimatore we are to use. (see compareEstimators)

estimator=5 #5-trioml,6-wang,7-lynchli,8-lynchrd,9-ritland,10-quellergt,11-dyadml
estimatorName="trioml" # change as needed with above
#Also choose which cutoffs to use. These are dependant on density plots
fshsCut<-0.4
hsurCut<-0.17

2011 Larval Relatedness Plots

main="2011 Larval snps"
fileName=paste("./outData/",main," coancestoryOutput", sep="")
load(file = fileName)
require(igraph)
relData <- output$relatedness[, c(2, 3, estimator)]
hist(relData[, 3], main = "Histogram of estimator")

lrelData <- relData
colnames(lrelData)[1] <- "from"
colnames(lrelData)[2] <- "to"
colnames(lrelData)[3] <- "weight"
lrelData$type <- "probRel"
relDataNoRows <- nrow(relData)
nrelData <- data.frame(matrix(ncol = 4, nrow = relDataNoRows))
colnames(nrelData) <- c("id", "name", "type", "label")
nrelData[, 1] <- relData[, 2]
nrelData[, 2] <- relData[, 2]
newrow = c(relData[1, 1], relData[1, 1], NA, NA)
nrelData = rbind(nrelData, newrow)
nrow(nrelData)
## [1] 1892
length(unique(nrelData$id))
## [1] 62
nrow(lrelData)
## [1] 1891
nrow(unique(lrelData[, c("from", "to")]))
## [1] 1891
nrelData <- data.frame(unique(nrelData[, 1:4]))
tst <- ifelse(grepl("^A", nrelData$id), nrelData$type <- "adult", nrelData$type <- "larvae")
nrelData$type <- tst
rm(tst)
larvNrel <- merge(nrelData, larv, by.x = "id", by.y = "LarvalRecords_LarvaID")
larvNrel <- larvNrel[, c(1:4, 66)]
head(nrelData)
##   id name   type label
## 1 71   71 larvae  <NA>
## 2 72   72 larvae  <NA>
## 3 73   73 larvae  <NA>
## 4 76   76 larvae  <NA>
## 5 75   75 larvae  <NA>
## 6 74   74 larvae  <NA>
head(larvNrel)
##    id name   type label YearOnly
## 1 100  100 larvae  <NA>     2011
## 2 101  101 larvae  <NA>     2011
## 3 103  103 larvae  <NA>     2011
## 4 104  104 larvae  <NA>     2011
## 5 105  105 larvae  <NA>     2011
## 6 107  107 larvae  <NA>     2011
net <- graph_from_data_frame(d = lrelData, vertices = larvNrel, directed = FALSE)
fileName = paste("./outData/", main, " net", sep = "", collapse = " ")
save(net, file = fileName)

Refine set and Plot relationships

fileName=paste("./outData/",main," net",sep="")
load(file = fileName)

#Limit to full siblings - this also makes it more sparse.
net.FS <- delete_edges(net, E(net)[weight<fshsCut])
l <- layout_with_fr(net.FS)

V(net.FS)$color=V(net.FS)$YearOnly #assign the "YearOnly" attribute as the vertex color
V(net.FS)$color=gsub("2011","indianred",V(net.FS)$color) #2011 will be red
V(net.FS)$color=gsub("2012","lightgoldenrod1",V(net.FS)$color) #2012 will be blue
V(net.FS)$color=gsub("2013","lightgreen",V(net.FS)$color) #2013 will be blue

E(net.FS)$weight<-E(net.FS)$weight*5

plot(net.FS, edge.arrow.size=0, edge.curved=0.2, vertex.size=5, vertex.color=V(net.FS)$color,vertex.frame.color="#555555",vertex.label=V(net)$name, vertex.label.color="black",vertex.label.cex=.7,edge.width=E(net.FS)$weight,layout=l)
#removed main=main;added edge.width=E(net.FS)$weight
title(paste(estimatorName,main),cex.main=3)
legend(x=-1.5, y=-1.1, c("2011","2012", "2013"), pch=21,col="#777777", cex=.8, bty="n", ncol=1)

2012 Larval Relatedness Plots

main="2012 Larval snps"
load(file = paste("./outData/",main," coancestoryOutput", sep=""))
require(igraph)
relData <- output$relatedness[, c(2, 3, estimator)]
hist(relData[, 3], main = "Histogram of estimator")

lrelData <- relData
colnames(lrelData)[1] <- "from"
colnames(lrelData)[2] <- "to"
colnames(lrelData)[3] <- "weight"
lrelData$type <- "probRel"
relDataNoRows <- nrow(relData)
nrelData <- data.frame(matrix(ncol = 4, nrow = relDataNoRows))
colnames(nrelData) <- c("id", "name", "type", "label")
nrelData[, 1] <- relData[, 2]
nrelData[, 2] <- relData[, 2]
newrow = c(relData[1, 1], relData[1, 1], NA, NA)
nrelData = rbind(nrelData, newrow)
nrow(nrelData)
## [1] 991
length(unique(nrelData$id))
## [1] 45
nrow(lrelData)
## [1] 990
nrow(unique(lrelData[, c("from", "to")]))
## [1] 990
nrelData <- data.frame(unique(nrelData[, 1:4]))
tst <- ifelse(grepl("^A", nrelData$id), nrelData$type <- "adult", nrelData$type <- "larvae")
nrelData$type <- tst
rm(tst)
larvNrel <- merge(nrelData, larv, by.x = "id", by.y = "LarvalRecords_LarvaID")
larvNrel <- larvNrel[, c(1:4, 66)]
head(nrelData)
##    id name   type label
## 1 151  151 larvae  <NA>
## 2 157  157 larvae  <NA>
## 3 156  156 larvae  <NA>
## 4 155  155 larvae  <NA>
## 5 154  154 larvae  <NA>
## 6 153  153 larvae  <NA>
head(larvNrel)
##    id name   type label YearOnly
## 1 134  134 larvae  <NA>     2012
## 2 135  135 larvae  <NA>     2012
## 3 136  136 larvae  <NA>     2012
## 4 137  137 larvae  <NA>     2012
## 5 138  138 larvae  <NA>     2012
## 6 139  139 larvae  <NA>     2012
net <- graph_from_data_frame(d = lrelData, vertices = larvNrel, directed = FALSE)
fileName = paste("./outData/", main, " net", sep = "", collapse = " ")
save(net, file = fileName)

Refine set and Plot relationships

fileName=paste0("outData/",main," net",sep=" ")
load(fileName) 

#Limit to full siblings - this also makes it more sparse.
net.FS <- delete_edges(net, E(net)[weight<fshsCut])
l <- layout_with_fr(net.FS)

V(net.FS)$color=V(net.FS)$YearOnly #assign the "YearOnly" attribute as the vertex color
V(net.FS)$color=gsub("2011","indianred",V(net.FS)$color) #2011 will be red
V(net.FS)$color=gsub("2012","lightgoldenrod1",V(net.FS)$color) #2012 will be blue
V(net.FS)$color=gsub("2013","lightgreen",V(net.FS)$color) #2013 will be blue

E(net.FS)$weight<-E(net.FS)$weight*5

plot(net.FS, edge.arrow.size=0, edge.curved=0.2, vertex.size=5, vertex.color=V(net.FS)$color,vertex.frame.color="#555555",vertex.label=V(net)$name, vertex.label.color="black",vertex.label.cex=.7,edge.width=E(net.FS)$weight,layout=l)
#removed main=main;added edge.width=E(net.FS)$weight
title(paste(estimatorName,main),cex.main=3)
legend(x=-1.5, y=-1.1, c("2011","2012", "2013"), pch=21,col="#777777", cex=.8, bty="n", ncol=1)

2013 Larval Relatedness Plots

main="2013 Larval snps"
load(file = paste("./outData/",main," coancestoryOutput", sep=""))
require(igraph)
relData <- output$relatedness[, c(2, 3, estimator)]
hist(relData[, 3], main = "Histogram of estimator")

lrelData <- relData
colnames(lrelData)[1] <- "from"
colnames(lrelData)[2] <- "to"
colnames(lrelData)[3] <- "weight"
lrelData$type <- "probRel"
relDataNoRows <- nrow(relData)
nrelData <- data.frame(matrix(ncol = 4, nrow = relDataNoRows))
colnames(nrelData) <- c("id", "name", "type", "label")
nrelData[, 1] <- relData[, 2]
nrelData[, 2] <- relData[, 2]
newrow = c(relData[1, 1], relData[1, 1], NA, NA)
nrelData = rbind(nrelData, newrow)
nrow(nrelData)
## [1] 8516
length(unique(nrelData$id))
## [1] 131
nrow(lrelData)
## [1] 8515
nrow(unique(lrelData[, c("from", "to")]))
## [1] 8515
nrelData <- data.frame(unique(nrelData[, 1:4]))
tst <- ifelse(grepl("^A", nrelData$id), nrelData$type <- "adult", nrelData$type <- "larvae")
nrelData$type <- tst
rm(tst)
larvNrel <- merge(nrelData, larv, by.x = "id", by.y = "LarvalRecords_LarvaID")
larvNrel <- larvNrel[, c(1:4, 66)]
head(nrelData)
##    id name   type label
## 1 186  186 larvae  <NA>
## 2 187  187 larvae  <NA>
## 3 188  188 larvae  <NA>
## 4 189  189 larvae  <NA>
## 5 190  190 larvae  <NA>
## 6 191  191 larvae  <NA>
head(larvNrel)
##    id name   type label YearOnly
## 1 185  185 larvae  <NA>     2013
## 2 186  186 larvae  <NA>     2013
## 3 187  187 larvae  <NA>     2013
## 4 188  188 larvae  <NA>     2013
## 5 189  189 larvae  <NA>     2013
## 6 190  190 larvae  <NA>     2013
net <- graph_from_data_frame(d = lrelData, vertices = larvNrel, directed = FALSE)
fileName = paste("./outData/", main, " net", sep = "", collapse = " ")
save(net, file = fileName)

Refine set and Plot relationships

fileName=paste("./outData/",main," net",sep="")
load(file = fileName)
#Limit to full siblings - this also makes it more sparse.
net.FS <- delete_edges(net, E(net)[weight<fshsCut])
l <- layout_with_fr(net.FS)

V(net.FS)$color=V(net.FS)$YearOnly #assign the "YearOnly" attribute as the vertex color
V(net.FS)$color=gsub("2011","indianred",V(net.FS)$color) #2011 will be red
V(net.FS)$color=gsub("2012","lightgoldenrod1",V(net.FS)$color) #2012 will be blue
V(net.FS)$color=gsub("2013","lightgreen",V(net.FS)$color) #2013 will be blue

E(net.FS)$weight<-E(net.FS)$weight*5

plot(net.FS, edge.arrow.size=0, edge.curved=0.2, vertex.size=5, vertex.color=V(net.FS)$color,vertex.frame.color="#555555",vertex.label=V(net)$name, vertex.label.color="black",vertex.label.cex=.7,edge.width=E(net.FS)$weight,layout=l)
#removed main=main;added edge.width=E(net.FS)$weight
title(paste(estimatorName,main),cex.main=3)
legend(x=-1.5, y=-1.1, c("2011","2012", "2013"), pch=21,col="#777777", cex=.8, bty="n", ncol=1)

2011-2013 Larval Relatedness Plots

main="2011-2013 Larval snps"
load(file = paste("./outData/",main," coancestoryOutput", sep=""))
require(igraph)
relData <- output$relatedness[, c(2, 3, estimator)]
hist(relData[, 3], main = "Histogram of estimator")

lrelData <- relData
colnames(lrelData)[1] <- "from"
colnames(lrelData)[2] <- "to"
colnames(lrelData)[3] <- "weight"
lrelData$type <- "probRel"
relDataNoRows <- nrow(relData)
nrelData <- data.frame(matrix(ncol = 4, nrow = relDataNoRows))
colnames(nrelData) <- c("id", "name", "type", "label")
nrelData[, 1] <- relData[, 2]
nrelData[, 2] <- relData[, 2]
newrow = c(relData[1, 1], relData[1, 1], NA, NA)
nrelData = rbind(nrelData, newrow)
nrow(nrelData)
## [1] 28204
length(unique(nrelData$id))
## [1] 238
nrow(lrelData)
## [1] 28203
nrow(unique(lrelData[, c("from", "to")]))
## [1] 28203
nrelData <- data.frame(unique(nrelData[, 1:4]))
tst <- ifelse(grepl("^A", nrelData$id), nrelData$type <- "adult", nrelData$type <- "larvae")
nrelData$type <- tst
rm(tst)
larvNrel <- merge(nrelData, larv, by.x = "id", by.y = "LarvalRecords_LarvaID")
larvNrel <- larvNrel[, c(1:4, 66)]
head(nrelData)
##    id name   type label
## 1 186  186 larvae  <NA>
## 2 187  187 larvae  <NA>
## 3 188  188 larvae  <NA>
## 4 189  189 larvae  <NA>
## 5 190  190 larvae  <NA>
## 6 191  191 larvae  <NA>
head(larvNrel)
##    id name   type label YearOnly
## 1 100  100 larvae  <NA>     2011
## 2 101  101 larvae  <NA>     2011
## 3 103  103 larvae  <NA>     2011
## 4 104  104 larvae  <NA>     2011
## 5 105  105 larvae  <NA>     2011
## 6 107  107 larvae  <NA>     2011
net <- graph_from_data_frame(d = lrelData, vertices = larvNrel, directed = FALSE)
fileName = paste("./outData/", main, " net", sep = "", collapse = " ")
save(net, file = fileName)

Refine set and Plot relationships

fileName=paste("./outData/",main," net",sep="")
load(file = fileName)
#Limit to full siblings - this also makes it more sparse.
net.FS <- delete_edges(net, E(net)[weight<fshsCut])
l <- layout_with_fr(net.FS)

V(net.FS)$color=V(net.FS)$YearOnly #assign the "YearOnly" attribute as the vertex color
V(net.FS)$color=gsub("2011","indianred",V(net.FS)$color) #2011 will be red
V(net.FS)$color=gsub("2012","lightgoldenrod1",V(net.FS)$color) #2012 will be blue
V(net.FS)$color=gsub("2013","lightgreen",V(net.FS)$color) #2013 will be blue

E(net.FS)$weight<-E(net.FS)$weight*5

plot(net.FS, edge.arrow.size=0, edge.curved=0.2, vertex.size=5, vertex.color=V(net.FS)$color,vertex.frame.color="#555555",vertex.label=V(net)$name, vertex.label.color="black",vertex.label.cex=.7,edge.width=E(net.FS)$weight,layout=l)
#removed main=main;added edge.width=E(net.FS)$weight
title(paste(estimatorName,main),cex.main=3)
legend(x=-1.5, y=-1.1, c("2011","2012", "2013"), pch=21,col="#777777", cex=.8, bty="n", ncol=1)

2011-2013 Larvae Relatedness with Adults

main="2011-2013 Larval snps with Adults"
load(file = paste("./outData/",main," coancestoryOutput", sep=""))
require(igraph)
relData <- output$relatedness[, c(2, 3, estimator)]
hist(relData[, 3])

lrelData <- relData
colnames(lrelData)[1] <- "from"
colnames(lrelData)[2] <- "to"
colnames(lrelData)[3] <- "weight"
lrelData$type <- "probRel"
relDataNoRows <- nrow(relData)
nrelData <- data.frame(matrix(ncol = 4, nrow = relDataNoRows))
colnames(nrelData) <- c("id", "name", "type", "label")
nrelData[, 1] <- relData[, 2]
nrelData[, 2] <- relData[, 2]
newrow = c(relData[1, 1], relData[1, 1], NA, NA)
nrelData = rbind(nrelData, newrow)
nrow(nrelData)
## [1] 34454
length(unique(nrelData$id))
## [1] 263
nrow(lrelData)
## [1] 34453
nrow(unique(lrelData[, c("from", "to")]))
## [1] 34453
nrelData <- data.frame(unique(nrelData[, 1:4]))
tst <- ifelse(grepl("^A", nrelData$id), nrelData$type <- "adult", nrelData$type <- "larvae")
nrelData$type <- tst
rm(tst)
larvNrel <- merge(nrelData, qslMetaLarvAndAdultsUnion, by = "id")
head(nrelData)
##    id name   type label
## 1 186  186 larvae  <NA>
## 2 187  187 larvae  <NA>
## 3 188  188 larvae  <NA>
## 4 189  189 larvae  <NA>
## 5 190  190 larvae  <NA>
## 6 191  191 larvae  <NA>
head(larvNrel)
##    id name   type label       pop YearOnly hatchDoY    lat    lon Site.Name catchDoY estimatedAge mother  father
## 1 100  100 larvae  <NA> M. peelii     2011 309.1240 -35.43 149.07 Murramore      326         9.39    MH6     FH6
## 2 101  101 larvae  <NA> M. peelii     2011 309.8434 -35.43 149.07 Murramore      326         8.67               
## 3 103  103 larvae  <NA> M. peelii     2011 297.4343 -35.43 149.07 Murramore      326        21.08               
## 4 104  104 larvae  <NA> M. peelii     2011 308.7643 -35.43 149.07 Murramore      326         9.75    MH6     FH6
## 5 105  105 larvae  <NA> M. peelii     2011 318.9451 -35.43 149.07 Murramore      340        11.73  Gilly Samwell
## 6 107  107 larvae  <NA> M. peelii     2011 309.6635 -35.43 149.07 Murramore      326         8.85    MH6     FH6
net <- graph_from_data_frame(d = lrelData, vertices = nrelData, directed = FALSE)
fileName = paste("./outData/", main, " net", sep = "", collapse = " ")
save(net, file = fileName)

Refine set and Plot relationships

fileName=paste("./outData/","2011-2013 Larval snps with Adults"," net",sep="")
load(file = fileName)
#Limit to full siblings - this also makes it more sparse.
net.FS <- delete_edges(net, E(net)[weight<fshsCut])
l <- layout_with_fr(net.FS)

plot(net.FS, edge.arrow.size=0, edge.curved=0, vertex.size=5,vertex.color=c("orange", "cyan")[(V(net.FS)$type=="larvae")+1], vertex.frame.color="#555555",vertex.label=V(net)$name, vertex.label.color="black",vertex.label.cex=.7,layout=l)
title(paste(estimatorName,main),cex.main=3)

2011-2013 Larvae with Parents

#Prepare data
prelDataA<-cbind(larvNrel[,c(1,13)],rep("1",nrow(larvNrel)),rep("mother",nrow(larvNrel)))

colnames(prelDataA)[1]<-"from"
colnames(prelDataA)[2]<-"to"
colnames(prelDataA)[3]<-"weight"
colnames(prelDataA)[4]<-"type"

prelDataB<-cbind(larvNrel[,c(1,14)],rep("1",nrow(larvNrel)),rep("father",nrow(larvNrel)))
colnames(prelDataB)[1]<-"from"
colnames(prelDataB)[2]<-"to"
colnames(prelDataB)[3]<-"weight"
colnames(prelDataB)[4]<-"type"
prelData<-rbind(prelDataA,prelDataB)
prelData[prelData==""] <- NA
rm(prelDataB);rm(prelDataA)
tmp<-prelData[complete.cases(prelData),]

#Create Vertices Data Frame
prelVert <- data.frame(matrix(ncol = 4, nrow = nrow(tmp)))
colnames(prelVert) <- c("id", "name", "type", "label")

prelVert$id<-tmp$to
prelVert$name<-tmp$to
prelVert$type<-tmp$type
prelVert<-rbind(prelVert, nrelData)
prelVert<-unique(prelVert)


prelData<-prelData[complete.cases(prelData),]

prelData <- subset(prelData, !from == 177) #these line needed to remove a few potential contaminants
prelData <- subset(prelData, !from == 314) #although I doubt there was actually contam as the otherhs are the same.
prelData <- subset(prelData, !from == 366)


require(igraph)
net.parents<-graph_from_data_frame(d=prelData, vertices=prelVert, directed=FALSE)

l <- layout_with_fr(net.parents)

#Colour vertices :parents and larvae
V(net.parents)$color=V(net.parents)$type #assign the "type" attribute as the vertex color then assign a colour based on that type.
V(net.parents)$color=gsub("mother","lightpink",V(net.parents)$color)
V(net.parents)$color=gsub("father","lightblue",V(net.parents)$color)
V(net.parents)$color=gsub("larvae","lightgreen",V(net.parents)$color)
V(net.parents)$color=gsub("adult","orange",V(net.parents)$color)

plot(net.parents, edge.arrow.size=0, edge.curved=0, vertex.size=5, vertex.color=V(net.parents)$color,vertex.frame.color="#555555",vertex.label=V(net.parents)$name, vertex.label.color="black",vertex.label.cex=.7,layout=l, main = "Full Siblings With Parents")

Session Info

all_labels()
##  [1] "Project_Template_and_Knitr"                   "Set_Global_Options"                           "unnamed-chunk-1"                              "unnamed-chunk-2"                             
##  [5] "Nominate Estimator to Use"                    "unnamed-chunk-3"                              "Larvae2011Net"                                "unnamed-chunk-4"                             
##  [9] "unnamed-chunk-5"                              "Larvae2012Net"                                "unnamed-chunk-6"                              "unnamed-chunk-7"                             
## [13] "Larvae2013Net"                                "unnamed-chunk-8"                              "unnamed-chunk-9"                              "Larvae2011-2013Net"                          
## [17] "unnamed-chunk-10"                             "unnamed-chunk-11"                             "unnamed-chunk-12"                             "LarvaeAdults2011-2013Net"                    
## [21] "unnamed-chunk-13"                             "unnamed-chunk-14"                             "unnamed-chunk-15"                             "HalfSibsColouredByYear"                      
## [25] "Include_Chunk_Labels_and_Session Information" "createNet"                                    "createNetA"
proc.time()-ptm
##    user  system elapsed 
##   17.21    2.37   33.97
#Session Information
sessionInfo()
## R version 3.3.0 (2016-05-03)
## Platform: i386-w64-mingw32/i386 (32-bit)
## Running under: Windows 7 x64 (build 7601) Service Pack 1
## 
## locale:
## [1] LC_COLLATE=English_Australia.1252  LC_CTYPE=English_Australia.1252    LC_MONETARY=English_Australia.1252 LC_NUMERIC=C                       LC_TIME=English_Australia.1252    
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
## [1] ggplot2_2.1.0       igraph_1.0.1        dplyr_0.4.3         ProjectTemplate_0.6 knitr_1.13         
## 
## loaded via a namespace (and not attached):
##  [1] Rcpp_0.12.5      digest_0.6.9     assertthat_0.1   plyr_1.8.4       grid_3.3.0       R6_2.1.2         gtable_0.2.0     DBI_0.4-1        formatR_1.4      magrittr_1.5     scales_0.4.0    
## [12] evaluate_0.9     stringi_1.1.1    lazyeval_0.1.10  rmarkdown_0.9.6  tools_3.3.0      stringr_1.0.0    munsell_0.4.3    yaml_2.1.13      parallel_3.3.0   colorspace_1.2-6 htmltools_0.3.5